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Abstract 

How best to design and redesign high-throughput experiments for zeohte synthe- 
sis is addressed. A model that relates materials function to chemical composition 
of the zeolite and the structure directing agent is introduced. Using this model, 
several Monte Carlo-like design protocols are evaluated. Multi-round protocols 
are found to be effective, and strategies that use a priori information about the 
structure-directing libraries are found to be the best. 
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1 Introduction 

Zeolites have found wide industrial applications in processes such as catalysis, molecular 
sieving, gas separation, and ion exchange. The significant catalytic activity and selectivity 
of zeolite materials are attributed to the large internal surface area and highly distributed 
active sites that are accessible through uniformly- sized pores. However, due to the desire to 
improve the quality of known zeolites that are small or faulted and the continuing quest for 
new zeolites with different functionalities, novel syntheses of zeolites with various frameworks 
and chemical compositions have drawn heavy research attention. To date, roughly 150 
framework structures have been reported. New structures are added at the pace of roughly 
6 per year. The use of organo-cation template molecules to provide structure direction has 
given rise to breakthroughs in zeolite science and provides much of the impetus for current 
synthetic efforts [|T], 0] . The addition of organic templates such as amines and alkylammonium 
ions to zeolite synthesis gels can affect the rate at which a particular material is formed and 
can make new structures or framework chemical compositions accessible. In some cases, a 
strong correlation between the geometry of the structure directing organic molecule and the 
zeolite pore architecture is found. The nature and extent of interactions between the organic 
templates and the inorganic components of a zeolite synthesis gel are key parameters that 
determine the final zeolite pore architecture 0, |]. 

We are interested in high-throughput synthesis of zeolites. High-throughput, or combi- 
natorial, methods allow for simultaneous creation of a large number of structurally diverse 
and complex compounds, generalizing the traditional techniques of single compound syn- 
thesis. Monte Carlo methods have been proposed and shown to be efficient methods for 
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library design and re-design in both material discovery ^] and small molecule design [0. 
Material discovery deals with continuous variables, such as composition variables and non- 
composition variables. Small molecule design deals with discrete variables, such as the iden- 
tities of template and ligand. For templated zeolite synthesis, we have both the continuous 
zeolite composition and non-composition variables and the discrete identity variables of the 
component parts of the organic structure-directing molecules. All of these variables affect 
the final zeolite material in an interrelated way. 

We here propose several strategies for templated, high-throughput zeolite synthesis. In 
section 2 we introduce a new random energy model that serves as a surrogate for experimental 
measurement of the figure of merit. This random energy model is based upon models used in 
material discovery and small molecule design. The random energy model is not fundamental 
to the protocols; it is introduced as a simple way to test, parameterize, and validate the 
various searching methods. In an experimental implementation, the random energy model 
will be replaced by the value returned by the screen. In section 3, we introduce several 
protocols for experimentation. In section 4, we compare results of each protocol and discuss 
some implications. Finally, we conclude in section 5. 

2 Random Energy Model 

To test proposed protocols in an efficient fashion, we need a model in lieu of the real 
experimental screening process. That is, we need a model that relates the chemical compo- 
sition of the zeolite and structure directing agent to the function, or figure of merit, of the 
material. This model is not essential to the protocols. It is simply a fast and cheap means 
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to test proposed protocols. This model will capture the essence of the physical system and 
provide validation for the protocols. Random energy models have proven useful for the study 
of several complex fields, including protein folding ^, protein molecular evolution [p!0| , 
material discovery and small molecule design 

To consider template-assisted zeolite synthesis, we use a model that combines ideas from 
both material discovery |^ and small molecule design 0]. The value of the figure of merit is 
naturally a sum of a zeolite energy and a well-chosen organic molecular energy: 



The zeolite framework composition variables Xi, d in number, are certainly key variables 
in i^zeoiite- The composition variables are specified by the mole fractions Xi, with < Xi < 1, 
and Y,i=iXi = 1- Typical elemental compositions in zeolite frameworks include silicon, 
oxygen, aluminum, phosphorus, boron, and germanium, among others. As discussed in 
the Appendix, the allowed compositions are constrained to lie within a simplex in d — 1 
dimensions. Several phases will exist for different compositions of the material. The figure 
of merit will be dramatically different between each of these distinct phases. To mimic this 
behavior, the composition variables are grouped in the Random Phase Volume Model into 
phases centered around points randomly placed within the allowed composition range. 
The phases form a Voronoi diagram, see Figure ^ The model is defined for any number of 
composition variables, and the number of phase points is defined by requiring the average 
spacing between phase points to be ^ = 0.25. To avoid edge effects, additional points are 
added in a belt of width 2,^ around the simplex of allowed compositions. The figure of merit 
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should change dramatically between composition phases. Moreover, within each phase a, 
the figure of merit should also vary with x — Xq, due to crystallinity effects such as crystallite 
size, intergrowths, defects, and faulting. Non-composition variables should also affect the 
measured figure of merit. The non-composition variables are denoted by the b variables Zj. 
Without loss of generality, each variable ranges —l<Zi<l. How the figure of merit changes 
with the non-composition variables should depend on the value of the composition variables. 
The non-composition variables fall within A^^ non-composition phases. There are a factor of 
10 fewer non-composition phases than composition phases. When x is in composition phase 
a and non-composition phase 7, the zeolite energy is given by 

^zeolite(x, z) = H - x", {G"}) + ^H{z, {G^}) (2) 

where {G"} and {G^} are phase-dependent random Gaussian variables with zero mean and 
unit variance. The {G"} are different in each composition phase, and the {G"' } are different 
in each non-composition phase. The form of the function H is 

H{w,, ...,wn, {G}) = Go + E E (3) 

k=l ii + ...-\-in=h 

n,...,i„>o 

The degree of the polynomial is chosen as 5 = 6. The set of coefficients {G} includes Gq and 
G^n...i„- The constant symmetry factors fi^...i^;k are given by 



The scale factor is chosen so that each term in the multinomial contributes roughly the 
same amount in the root-mean-square sense. For the composition variables C,h = ^/2, and 
for the non-composition variables = / {z'^)Y^^ = (3/7)^/^. The an are chosen so that 
the multinomial terms contribute 40% as much as the corresponding constant, phase terms, 
Go in eq ^ in the root-mean-square sense. 

The organic structure-directing molecules are uniquely characterized by composition, 
such as the identity of the template and ligands, as shown in Figure Note that template 
here denotes the framework of the structure-directing molecule, onto which ligands are at- 
tached, not the entire structure-directing molecule. For simplicity, we consider the small 
molecule to consist of one template and six binding ligands, each drawn from a fragment li- 
brary. To quantitatively describe each template and ligand, six weakly correlated descriptors 
are used. Their numeric values specify the characteristics of each template or ligand. 

We again use a random energy model for the molecules to account for the contributions 
to structure direction arising from interactions between the zeolite and template and from 
interactions between the zeolite and each of the ligands. In addition, synergistic effects be- 
tween the ligands and template are incorporated. Consider, for example, a molecule made 
from template number m from the template library and ligands number Si, . . . ,Sq from the 
ligand library. The template is characterized by six descriptors, Di^\ . . . ,Dq"\ Similarly, 
each ligand is characterized by six descriptors, di^\ . . . , dQ^\ We denote the template con- 
tribution to binding by E^^ and the ligand contributions by Ei^. We denote the contribution 
due to synergistic ligand-ligand interactions by -Ell and the contribution due to synergistic 
template-ligand interactions by -Etl- The total molecular contribution to the figure of merit 
is, then, 
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-E'molecule = -^l + -Et + -E'lL + -E'tL (5) 

Each of these terms is given in the form of a random polynomiah 



£;l(x,z) = ^F[4^\...,4'^\{Gl(x,z)}], Gl(x,z) = A.,iG^ + ^A,,iGZ (6) 

i=l 

E^{^,z) = F[Dt\...,Dt\{Gr{^,z)}], Gt(x,z) = A.,2G'^ + ^A,,2G^ (7) 

Ell(x, z) = ^ /^.(x, z)F[4';\ 4';\ 4';), 4f 4f 4f {Gll(x, z)}], 

i=l 

GLL(x,z) = A,,3Gi:^ + iA,,3Gi:Z (8) 
i?TL(x,z) = j:h,i^,z)F[4l\4;\4;\Dt:\Dl:\D^-\{G^^i^,z)}l 

i=l 

Gtl(x, z) = X,^,G'il + ^K,,Gi^l (9) 



where individual descriptors from each template and ligand are given in corresponding li- 
braries. The building block for our random energy model as a function of those descriptors 



IS 



F{w,, {G}) = E E . . . < (10) 



fc=0 ii + ... + in = fc 
il,...,in>0 



The {G}, again, will be composed of a composition and a non-composition piece, as indicated 
in eqs We use g = 6 and n = 6. Since the values of the descriptors, Wi, will be drawn 
from a Gaussian random distribution of zero mean and unit variance, we set the scaling 
factor, ^F, by 

1 1 

9-2 

(11) 



In general the polynomial coefficients, {G} in eq (^^, are functions of both zeolite variables 



X and z. However, for simplicity, we assume that the coefficients are phase-dependent only. 
Both types of zeolite phases have a significant impact on the energy function. The variables 
{G"}, {Gl}, {G^}, {Gt}' {^ll}' {^ll}' {^tl}' {G^l} are fixed random Gaussian vari- 
ables with zero mean and unit variance. These coefficients are different in each phase. The 
structural constant hi indicates the strength of the interaction at structure directing site i. 
This interaction function depends on the zeolite variables x and z as 

/^.(x, z) = - x", {Cn) + Ih{z, {G]}) (12) 

As before {G"} and {G]} are sets of phase-dependent fixed random Gaussian variables 
with zero mean and unit variance. The parameter an is adjusted so that the multinomial 
terms contribute roughly the same as the constant, phase terms on average. In -Ell? only 
synergistic interactions between neighboring ligands are considered, and it is understood 
that l-r refers to h in eq |^. In principle, the polynomial in eq ^ could be a function of all 
12 descriptors of both ligands. We assume, however, that important contributions come 
from interactions among three randomly chosen distinct descriptors of ligand li, dj'^\ d^''^\ 
and d^^^\ and another three randomly chosen distinct descriptors of ligand Zj+i, dj^^^\ dj^^^\ 
and c^jg^^^- Similarly, we assume that template-ligand contributions come from interactions 
between three randomly chosen distinct descriptors of the ligand, d^^l\df^\ and d^j^!^\ and 
another three randomly chosen distinct descriptors of the template, D^^\d^^\ and D^^\ 
Both ji and ki are descriptor indices ranging from 1 to 6. We assume that these indices 
depend only on the template, since the molecular details of the zeolite have been suppressed. 
The constants \x,i are adjusted so that the synergistic terms in compositional phases will 

8 



contribute on average as Ei^ : E'y : i^LL : -E'tl = 1 : 1 : 2 : 1.2. Furthermore, we adjust 
the Xx^i so that the total molecular contribution from the interaction with the compositional 
variables accounts for roughly 15% of the constant, compositional phase term in -Ezeoiite in the 
root-mean-square sense. We similarly adjust the values of X^^i in non-compositional phases 
so that the contribution from the interaction with the non-composition variables is also 15%. 

Finally, we minimize the energy returned by the model. That is, we sample the figure 
of merit, E, by exp {—[3E). 

3 Protocols for high-throughput experimentation 

As in the small molecule case 0, we first build the template and ligand libraries. We 
denote the size of the template library by A^t = 15 and the size of the ligand library by 
Al = 1000. In a real experiment, the six descriptors would then be calculated for each 
template and ligand. In the simulated experiment, the values of the six descriptors of each 
ligand and template are extracted from a Gaussian random distribution with zero mean and 
unit variance. In the simulated experiment, we also associate two sets of random interaction 
descriptor indices to each template for the interaction terms in eq ^ and |^. 

For comparison, we fix the total number of samples to be synthesized at 100000 for all 
protocols. That is, we keep roughly the same experimental cost. 

We consider several single pass protocols. Current high-throughput experiments uni- 
formly tend to perform a single-pass, grid search on all continuous variables. To mimic 
current experiments, we use such a grid search on the composition and non- composition 
variables. For the discrete molecular variables, we instead randomly pick one template and 
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six ligands from the corresponding libraries. In another single-pass protocol, which we call 
random, we search randomly over all variables, i. e. we choose the composition variables and 
non-composition variables at random as well. 

We also consider several multi-pass protocols. Unlike single-pass protocols, multi-pass 
protocols allow us to learn about the system as the experiment proceeds. In analogy to Monte 
Carlo methods, each round of combinatorial chemistry experiments corresponds to a move 
in a Monte Carlo simulation. Instead of tracking one system with many configurational 
degrees of freedom, however, many samples are tracked. This results in a rather diverse 
population of samples and increases the opportunity for a few zeolite samples to survive 
more elaborate tests for application performance, tests that are only roughly correlated with 
the primary screen. The existence of these secondary and primary screens is the reason why 
sampling the figure of merit to find several promising compounds is so important. We use a 
Monte Carlo protocol so as to maintain diversity in the samples. Genetic algorithms, while 
powerful at local optimization, reduce the diversity of the population in the cross-over step 
and tend to fail diversity tests 0. To be consistent with current experiments, we prefer to 
synthesize 1000 samples simultaneously and perform Monte Carlo moves on every sample 
for 100 rounds. The initial sample configurations are assigned by the random protocol. 

We change the zeolite variables as well as the components of the structure directing 
agent at each round of the Monte Carlo protocol. In the simplest approach, we perturb both 
the composition variables and the non-composition variables of the zeolite part around the 
original value using the traditional Metropolis-type method. For the non-composition vari- 
ables, z, periodic boundary conditions are used. For the composition variables, x, reflecting 
boundary conditions are used, as discussed in the Appendix. For the molecular part, either 
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the template is changed with probabihty ^template, or one of the six hgands is picked randomly 
to change in each Monte Carlo move. The new configurations are updated according to the 
acceptance rule at /3, the inverse of the protocol temperature. All the samples are sequen- 
tially updated in one Monte Carlo round. For the molecular part, either a simple Metropolis 
method or a biased Monte Carlo method can be incorporated. In the simple Metropolis 
method, the new template or ligand is selected at random with a uniform probability from 
the corresponding libraries. The final acceptance rule for an unbiased protocol, which is 
applied after the zeolite and molecular variables have been modified, is 

acc(o n) = min[l, exp{—pAE)] (13) 

The biased Monte Carlo schemes allow us to generate a trial configuration with a prob- 
ability that depends on the potential energy of that configuration. These schemes have been 
proven to be highly successful in small molecule high-throughput experiment design 
where biased energy forms are constructed from either theory or pre-experiments. In the 
present case, since the total molecular energy is a function of both the composition and 
non-composition variables, it is not feasible to construct the bias from pre-experiments. We, 
therefore, use the theoretical bias directly from the random energy model. For a zeolite at 
compositional phase a and non-compositional phase 7, the bias energy form for template m 
is 

E(™)(x, z) = F[D^r\ ...,Dt\ {Gt(x, z)}], Gt(x, z) = A^.^G? + ^A^.^G^ (14) 
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And the bias for hgand i is 



e«(x, z) = F[d?, {Gl(x, z)}], GlIx, z) = A^.^G^ + ^A^.^GZ 



(15) 



These biases tend to exhibit a large gap between a few dominant templates and ligands 
and the rest. To ensure the participation of more ligands and templates in the strategy, we 
introduce cutoff energies for the hgand and template, Cc and Ec- We choose Cc to be the 21st 
lowest ligand energy and Ef. to be the 4th lowest template energy. With this cutoff energy, 
the biased energy, ■* , for the i*'* ligand becomes 



In one Monte Carlo move, if the i*" ligand in the molecule is to be changed at new 
zeolite values x' and z', the biased probability for selecting another ligand from the ligand 
library is 




(16) 



And the biased energy, E^' , for the m*'* template become 




(17) 



(18) 



E;=\exp[-/3e[,^^(x',z')] 
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For reverse move, we have 



f[E{o)] = -M-pei^'M] (19) 
E&exp[-/?e(^)(x,z)] 



Similarly, if the template is desired to be changed at new zeolite values x' and z', the 
biased probability for selecting template m' is 



E,-Jiexp[-/5E^''(x',z')] 



For reverse move, we have 



exp[-(3E!r\^,z)] 
EfJiexp[-/54^'Hx,z)] 



/[^(°)] = ' r .X (21) 



To satisfy the detail balance, the acceptance rule becomes 



acc(o^n) =min|l,|||Mexp(-/3AE)| (22) 



Parallel tempering is known to be a powerful tool for searching rugged energy land- 



scapes [11, 12, 13]. Parallel tempering can be combined with both Metropolis Monte Carlo 
and biased Monte Carlo. In parallel tempering, the samples are divided into k groups. 
The first group of samples is simulated at the second group is at P2, and so on, with 
Pi < P2 < ■ ■ ■ < Pk- At the end of each round, samples in group i < k are allowed to 
exchange configurations with samples in group i + 1 with probability Pex- The acceptance 
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rule for a parallel tempering exchange move is 

acc(o n) = min[l, exp(A/3 AE)] (23) 

where A(5 = jSi — /Sj+i and AE is the difference in energy between the sample in group i 
and the sample in group i + 1 before the exchange is made. It is important to notice that 
this exchange step is experimentally cost-free. Nonetheless, this step can be dramatically 
effective at facilitating the protocol to escape from local minimum. The number of groups, 
the number of samples in each group, the value of /Jj, and the exchange probability, Pex, are 
experimental parameters to be tuned. 

4 Results and Discussion 

The size of the library is fixed at A'^r = 15 and A^l = 1000. We adjust Xx^i and \z,i so that 
the synergistic terms will contribute on average as Ei^ : E^^ : : -S'TL = 1 : 1 : 2 : 1.2 in the 
composition and non-composition phases. We also want a model in which the contribution 
from zeolite and molecule will be roughly of the same order for the optimal configurations 
obtained by grid or random protocols. Due to the complicated molecular terms, we find the 
parameters so that this occurs by trial and error. We first adjust all the au values, then 
we adjust the \x,i and \z,i in eqs As mentioned, we adjust the \x,i and \z,i so that 
the total molecular contribution from either type of phase accounts for roughly 15% of the 
corresponding constant phase term in i^zeoUte- 

To locate optimal parameters for the protocols, we perform a few short pre-experiments. 
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Since the size of the template hbrary is relatively small, the optimal value for the probability 
of changing a template is ^template = 0.02. The maximum random displacements are |Ax| = 
0.1/\/d — 1 and | Az| = 0.2 in the composition space and non-composition space, respectively. 
For simple Metropolis Monte Carlo, the optimal inverse protocol temperature is /3 = 50 for 
X3 and X4 and f3 = 20 for X5. For the biased Monte Carlo schemes, the optimal inverse 
protocol temperature is (3b = 500 for X3 and X4 and (3b = 200 for X5. For biased Monte 
Carlo combined with parallel tempering, it is optimal to have the samples divided into three 
subsets, with 25% of the population at (3i = ^(3b, 50% at P2 = Pb, and 25% at P3 = 2(3b. 
The switching probability, Pex? is 0.1. Determination of these parameter values corresponds 
experimentally to gaining familiarity with the protocol on a new system. 

All the protocols are tested with increasing numbers of composition and non- composition 
variables. Results are shown in fig ^. From this figure, one can see that all multi-round 
Monte Carlo protocols are better than the single-pass protocols such as grid and random. 
The simple Metropolis method can find optimum samples that are twice as superior as those 
from the grid or random protocols. Most importantly, the biased Monte Carlo is especially 
efficient. The optimum sample energies from biased Monte Carlo far exceed the optimum 
sample energies from simple Metropolis with even 1000 rounds. Biased Monte Carlo is able to 
find such favorable samples mainly by finding much more favorable designs for the structure- 
directing agent. Parallel tempering is noticeably effective only for the more complicated X4 
and X5 systems. For the relatively simple systems, such as X3, it is more important to keep 
all 1000 of the samples at the optimum temperature. 
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5 Conclusion 

We have presented a model for high-throughput zeohte synthesis. As in previous studies, 
multi-pass Monte Carlo methods work better than do single-pass protocols. Sophisticated 
biased Monte Carlo schemes are highly efficient and much better than simple Metropolis 
Monte Carlo. Parallel tempering is the best method for complicated systems with 5 or more 
framework chemical compositional variables. 
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Appendix: reflecting boundary conditions 

The valid region of composition variables form a simplex, since they are non-negative 
and sum to unity. For those points near the corner of the allowed simplex, a high percentage 
of proposed moves will go beyond the allowed region. If those moves are simply rejected, 
the composition point will get stuck in the corner for a long time, and the limited sampling 
chances will be wasted. We want an algorithm that satisfies detail balance while avoiding 
getting stuck in the corners. 

The scheme we have chosen ensures that all moves from any point in the simplex lead 
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to new valid points in the simplex. The scheme basically involves many reflections over the 
d — 2 dimensional hyper-planes that define the boundaries of the simplex. Figure ^ shows 
how this works in the case of ci = 3 for a single refiection. We first define the unit normals, 
rij, to each of the faces of the simplex. We also find the constants Cj that define the faces by 
the equation x ■ rij = q. In the case of figure ^, there are three such unit normals and three 
such constants. We initially define y = x(o). The algorithm consists of the following steps: 

1. Test whether the proposed new point, x, is in the allowed region. If so, define the new 
point x(n) to be x and stop. 

2. If not, find the face i for which the quantity = (q — y ■ ni)/[(x — y) ■ iij] is minimal, 
taking into account only those faces for which both ti and (x — y) ■ rij are positive. 

3. Define y ^ y + tj(x - y). 

4. Refiect the point through this face by the equation x ^ x — 2(x ■ n,; — Cj)nj. 

5. Continue with step 1. 

This algorithm converges, because it decreases the magnitude of x at each step, until x is 
within the allowed region. The algorithm is also reversible. That is, for each forward move 
x(o) — > X that maps to the new position x(n), there is always a reverse move x(n) — x' that 
maps to the original position x(o). Furthermore, when the original move is in the allowed 
move sphere, |x(o) — x| < Ax, the reverse move is as well, |x(n) — x'| < Ax. Given that we 
impose the detailed balance condition, then, these reflecting boundary conditions provide 
a valid Monte Carlo move from an old point in the allowed simplex to a new point in the 
allowed simplex. 
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These reflecting boundary conditions can be interpreted in a geometrical way. Essen- 
tially, the reflections perform a type of billiards in the d — 1 dimensional simplex. In each 
attempted move, we chose a displacement x— x(o). We can interpret the composition variable 
as a small ball that moves at a constant speed along this trajectory. The motion is contin- 
ued until either a boundary is encountered or the entire length of the trajectory has been 
traveled. If a boundary of the simplex is encountered, the ball reflects off the hyper-plane by 
Newtonian mechanics, and the trajectory continues along the new, reflected direction. The 
motion, including possible additional reflections, is continued until the ball has traveled a 
distance equal to the chosen length of the displacement, |x — x(o)|. The location of the ball 
at the end of the trajectory then gives the new point x(n). The reverse move is given by 
exactly the reverse motion, which has the same total displacement. 
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Figure 1: Schematic of the random energy model for templated zeohte synthesis. 
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Figure 2: The minimum figure of merit found by different protocols witli different number 
of composition (x) and non-composition (z) variables. Tlie complexity of tlie structure- 
directing agent libraries is the same in all cases. The results are scaled to the minimum 
found by the grid searching method. Each value is averaged over scaled results on 10 different 
instances. 




Figure 3: Schematic of the reflecting boundary conditions used to obtain valid composition 
variables. 
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